#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#define INF 10000

int main(int argc, char *argv[]) {
	float **matriz;
	float **matrizAEstrela;
	float *matrizB;
	float *matrizBEstrela;
	float *X0;
	float *auxMultiplicacao;
	float *X1;
	int i, j, iteracao = 1;
	float row_test = 0.0;

	int J_ORDER;
	int J_ROW_TEST;
	int J_ITE_MAX;
	float J_ERROR;
	
	struct timeval start, end;
	int sec, usec;
	
	/*
	 * Abertura do arquivo
	 */
	FILE *fp = fopen(argv[1], "r");
	if (fp == NULL) {
		printf("Erro: abertura do arquivo!\n");
		exit(EXIT_FAILURE);
	}
	rewind(fp);
	/*
	 * Leitura da ordem da matriz, da linha de teste, do erro e
	 * do numero maximo de iteracoes
	 */
	fscanf(fp, "%d%d%f%d", &J_ORDER, &J_ROW_TEST, &J_ERROR, &J_ITE_MAX);

	
	/*
	 * Alocacao das matrizes e vetores
	 */
	matriz = (float **) malloc(J_ORDER * sizeof(float *));
	matrizAEstrela = (float **) malloc(J_ORDER * sizeof(float *));
	matrizB = (float *) malloc(J_ORDER * sizeof(float));
	matrizBEstrela = (float *) malloc(J_ORDER * sizeof(float));
	X0 = (float *) malloc(J_ORDER * sizeof(float));
	auxMultiplicacao = (float *) malloc(J_ORDER * sizeof(float));
	X1 = (float *) malloc(J_ORDER * sizeof(float));
	
	for (i = 0; i < J_ORDER; i++) {
		matriz[i] = (float *) malloc(J_ORDER * sizeof(float));
		matrizAEstrela[i] = (float *) malloc(J_ORDER * sizeof(float));
	}
	
	gettimeofday(&start, NULL);

	/*
	 * Leitura da matriz a partir do arquivo
	 */
	printf("Reading matrices... ");
	for (i = 0; i < J_ORDER; i++)
		for (j = 0; j < J_ORDER; j++) {
			fscanf(fp, "%f", &matriz[i][j]);
			matrizAEstrela[i][j] = matriz[i][j];
		}

	for (i = 0; i < J_ORDER; i++) {
		fscanf(fp, "%f", &matrizB[i]);
		matrizBEstrela[i] = matrizB[i];
	}

	fclose(fp); /* Fecha arquivo */
	
	/*
	 * Divide cada linha pelo elemento da diagonal nas matrizes: 
	 * matrizBEstrela e matrizAEstrela.
	 * Na matrizAEstrela multiplica por -1
	 */
	printf("Done.\nCalculating... ");
	for (i = 0; i < J_ORDER; i++) {
		matrizBEstrela[i] /= matriz[i][i];
	}


	for (i = 0; i < J_ORDER; i++) {
		for (j = 0; j < J_ORDER; j++) {
			if (i != j) {
				matrizAEstrela[i][j] /= matriz[i][i];
				matrizAEstrela[i][j] *= -1;
			}
			matrizAEstrela[i][i] = 0;
		}
	}

	/*
	 * Calcula X(k+1) = matrizAEstrela * X(k)+ matrizBEstrela
	 * Ate que o numero maximo de iteracoes seja atingido ou 
	 * o erro estaja dentro do especificado
	 */
	while (iteracao < J_ITE_MAX && fabs(matrizB[J_ROW_TEST] - row_test) > J_ERROR) {

		for (i = 0; i < J_ORDER; i++)
			auxMultiplicacao[i] = 0;
		
		/* auxMultiplicacao = matrizAEstrela * X(k) */
		for (i = 0; i < J_ORDER; i++)
			for (j = 0; j < J_ORDER; j++)
				auxMultiplicacao[i] += matrizAEstrela[i][j] * X0[j];

		/* X1 = auxMultiplicacao + matrizBEstrela */
		for (i = 0; i < J_ORDER; i++)
			X1[i] = auxMultiplicacao[i] + matrizBEstrela[i];

		/* Calculo do Erro */
		row_test = 0;
		for (j = 0; j < J_ORDER; j++) {
			row_test += matriz[J_ROW_TEST][j] * X1[j];
		}

		// X0 = X1
		for (i = 0; i < J_ORDER; i++) {
			X0[i] = X1[i];
		}

		iteracao++;
	}
	
	gettimeofday(&end, NULL);
	sec = end.tv_sec - start.tv_sec;
	usec = fabs(end.tv_usec - start.tv_usec);
	printf("Done.\nEXECUTION TIME: %d.%d seconds\n\n", sec, usec);

	printf("---------------------------------------------------------\n");
	printf("Iterations: %d\n", iteracao);
	printf("RowTest: %d => [%f] =? %f\n", J_ROW_TEST, row_test, matrizB[J_ROW_TEST]);
	printf("---------------------------------------------------------\n");
	
	/* Libera espacos alocados para matrizes e vetores */
	for (i = 0; i < J_ORDER; i++) {
		free(matriz[i]);
		free(matrizAEstrela[i]);
	}
	free(matriz);
	free(matrizAEstrela);
	free(matrizB);
	free(matrizBEstrela);
	free(X0);
	free(auxMultiplicacao);
	free(X1);

	return 0;
}
